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Abstract 

This paper presents a methodology to determine the parameters used in the sim- 
ulation of delanrination in composite materials using decohesion finite elements. A 
closed-form expression is developed to define the stiffness of the cohesive layer. A 
novel procedure that allows the use of coarser meshes of decohesion elements in 
large-scale computations is proposed. The procedure ensures that the energy dissi- 
pated by the fracture process is correctly computed. It is shown that coarse-meshed 
models defined using the approach proposed here yield the same results as the 
models with finer meshes normally used in the simulation of fracture processes. 


1 Introduction 


Delamination, or interfacial cracking between composite layers, is one of the 
most common types of damage in laminated fibre-reinforced composites due 
to their relatively weak interlaminar strengths. Delamination may arise un- 
der various circumstances, such as low velocity impacts, or bearing loads in 
structural joints. The delamination failure mode is particularly important for 
the structural integrity of composite structures because it is difficult to detect 
during inspection. 

The simulation of delamination using the finite element method (FEM) is 
normally performed using the Virtual Crack Closure Technique (VCCT) [1], 
or decohesion finite elements [2]- [3]. 



The VCCT is based on the assumption that the energy released during de- 
lamination propagation equals the work required to close the crack back to its 
original position. Based on this assumption, the single-mode components of 
the energy release rate are computed from the nodal forces and nodal relative 
displacements [1]. Delamination growth is predicted when a combination of 
the components of the energy release rate is equal to a critical value. 

There are some difficulties when using the VCCT in the simulation of pro- 
gressive delamination. The calculation of fracture parameters requires nodal 
variable and topological information from the nodes ahead and behind the 
crack front. Such calculations are tedious to perform and may require remesh- 
ing for crack propagation. 

The use of decohesion finite elements can overcome some of the above difficul- 
ties. Decohesion finite elements can predict both the onset and non-self-similar 
propagation of delamination. However, the simulation of progressive delam- 
ination using decohesion elements poses numerical difficulties related with 
the proper definition of the stiffness of the cohesive layer, the requirement 
of extremely refined meshes, and the convergence difficulties associated with 
problems involving softening constitutive models. 

This work addresses the proper definition of interface stiffness and proposes a 
novel procedure to allow the use of coarse meshes in the simulation of delam- 
ination. 

A brief description of the kinematics and constitutive model of a previously 
proposed decohesion element [2]- [3] is presented. A closed-form expression is 
developed, replacing the empirical definitions of the stiffness of the cohesive 
layer that are normally used. A solution to use coarse meshes in the simulation 
of delamination is proposed. It is demonstrated that the proposed solution can 
yield results as accurate as the ones obtained using very refined meshes. 


2 Cohesive zone model approach 


The Cohesive Zone Model (CZM) approach [4]— [6] is one of the most commonly 
used tools to investigate interfacial fracture. The CZM approach assumes that 
a cohesive damage zone, or softening plasticity, develops near the crack tip. 

The CZM links the microstructural failure mechanism to the continuum fields 
governing bulk deformations. Thus, a CZM is characterized by the properties 
of the bulk material, the crack initiation condition, and the crack evolution 
function. 
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Cohesive damage zone models relate traction to displacement jumps at an in- 
terface where a crack may occur. Damage initiation is related to the interfacial 
strength, i.e. , the maximum traction on the traction-displacement jump rela- 
tion. When the area under the traction-displacement jump relation is equal to 
the fracture toughness, the traction is reduced to zero, and new crack surfaces 
are formed. 


2.1 Kinematics and constitutive relation of cohesive zone models 


The cohesive zone model used here was previously proposed by the authors 
[2]- [3]. A brief description of the model, with special focus on the kinematics 
and constitutive damage model, is presented. 

Consider a domain D containing a material discontinuity, T d , which divides 
the domain D into two parts, D+ and D-, as shown in Figure 1. 



Fig. 1. Body D crossed by a material discontinuity T d in the undeformed configu- 
ration. 

Prescribed tractions, C, are imposed on the boundary T F , whereas prescribed 
displacements are imposed on the boundary T u . The stress field inside the 
domain, cr ^ , is related to the external loading and the closing tractions r ? + , t~ 
in the material discontinuity through the equilibrium equations: 


< Tijj — 0 in (1) 

(TijUj = ti on T F (2) 

= r+ = —t~ = (Jijrif on T d (3) 

The displacement jump across the interface of the material discontinuity re- 
quired in the constitutive model, [u ?; ] , can be obtained as a function of the 
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displacement of the points located on the top and bottom side of the interface, 
uf and u~ respectively: 


[«*! = u t - u i ( 4 ) 

where uf are the displacements with respect to the fixed Cartesian coordinate 
system. A co-rotational formulation is used in order to express the compo- 
nents of the displacement jumps with respect to the deformed interface. The 
coordinates Xi of the deformed interface can be written as [7]: 


x { 




( 5 ) 


where X t are the coordinates of the undeformed interface. 

The components of the displacement jump tensor in the local coordinate sys- 
tem on the deformed interface, A m , are expressed in terms of the displacement 
field in global coordinates: 


A. m — [m ;J (6) 

where Q„ u is the rotation tensor. 

The constitutive operator of the interface, D ]t , relates the element tractions, 
Tj, to the displacement jumps, Ay 


T j ~ Djj Ai (7) 

A requirement of the constitutive relations of cohesive zone models is that the 
energy dissipated in the process of fracture needs to be computed accurately. 

Under single-mode loading, controlled energy dissipation is achieved assur- 
ing that the area under the traction- displacement jump relation equals the 
corresponding fracture toughness, as illustrated in Figure 2. 

Under mixed-mode loading, a criterion established in terms of an interaction 
between components of the energy release rate needs to be used. 

There are several types of constitutive operators presented in the literature, 
depending on the constitutive equations used for the simulation of the de- 
lamination: Tvergaard and Hutchinson [8] used a trapezoidal law, Cui and 
Wisnom [9] proposed a perfectly plastic rule, Needleman used a polynomial 
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Fig. 2. Constitutive equations under Mode I and Mode II loading. 


law, [10], and later an exponential law [11]. In this paper, a bilinear constitu- 
tive equation is used. 

The constitutive damage model used here, formulated in the context of Dam- 
age Mechanics (DM), was previously proposed by the authors [2], [3]. All the 
details of the constitutive model are presented in reference [2] and [3] and will 
not be repeated here. 

The constitutive model prevents interpenetration of the faces of the crack 
during closing, and a Fracture Mechanics-based criterion is used to predict 
crack propagation. The formulation of the damage model is summarized in 
Table 1, where and ijp are the free energy per unit surface of the damaged 
and undamaged interface, respectively. Sij is the Kronecker delta, and d is 
a scalar damage variable. The parameter A is the norm of the displacement 
jump tensor (also called equivalent displacement jump norm), and it is used to 
compare different stages of the displacement jump state so that it is possible to 
define such concepts as ‘loading’, ‘unloading’ and ‘reloading’. The equivalent 
displacement jump is a non-negative and continuous function, defined as: 


^ — \Z(A 3 ) 2 + (A s h ear ) 2 (8) 

where (•) is the MacAuley bracket defined as (x) = | (x + |x|). A 3 is the 
displacement jump in mode I, i.e., normal to midplane, and A s ^ ear is the 
Euclidean norm of the displacement jump in mode II and in mode III: 

A s h. e ar = \j (Al) 2 + (A 2 )" (9) 
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Table 1. Definition of the constitutive model. 


Free Energy 

V’(A.d) = (1 — d) ip° (A i 

) - dip 0 (S 3i {- 

-A3)) 

Constitutive equation 

Ti = Si = (! - d) SijKAj - d SijKhj 

(-As) 

Displacement jump norm 

X — \J (A3} 2 + {A s h ear ) 2 



Damage criterion 

F(X t ,r t ) := G (A 4 ) - G ( 

?’*) <0 Vt > 0 


C (\)- A/ ( A - A °) 
b ^ A ' ) “ A(A/-AO) 



Evolution law 

j _ ,-dF(\r) _ • 9G(A) . 

a ~ r d A — ^ d\ ’ 

r = pL 


Load/unload conditions 

l'i > 0 ; F (A*, r*) < 0 

; fiF ( A*, r-*) 

= 0 


r* = max { r° , max s A s } 

0 < 8 < t 



The evolution of damage is defined by G(-), a suitable monotonic scalar func- 
tion, ranging from 0 to 1. fi is a damage consistency parameter used to define 
loading-unloading conditions according to the Kuhn- Tucker relations. r l is the 
damage threshold for the current time, t, and r° denotes the initial damage 
threshold. 


A 0 is the onset displacement jump, and it is equal to the initial damage thresh- 
old r°. The initial damage threshold is obtained from the formulation of the 
initial damage surface or initial damage criterion. A? is the final displacement 
jump, and it is obtained from the formulation of the propagation surface or 
propagation criterion [2]- [3]. 


The formulation proposed allows an explicit integration of the constitutive 
model. The algorithm is implemented as shown in Table 2. 
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Table 2. Algorithm of the constitutive model. 


Initial data for time t+1 

Material properties: G IC , G n c, G inc , E , //, T s ° ftear , t 3 ° 
Current values: A 3 , A shear, d f 


1. Mixed-mode ratios: 


P = 


( A 3 )+A ahear ’ 


B = 


T 


l+2/3 2 -2 /3 


2. Pure mode onset displacement jump: 


A 0 = h- 

I< 


i = 3, shear 


3. Mixed-mode onset displacement jump: A 0 = J (A!]) + (A^ ear ) — (A^) [B] 


4. Mixed-mode final displacement jump: A^ = T [ Gj c + (Gjj c — G/ c ) [B] 


5. Displacement jump norm: A— \/ (^3) 2 + (A sh.earY 


6. Update internal variables: r f = AO i 


• r t+l = max 


{r*, A} 


i+1 A/(r«+*-A°) 

r t+i(Af-A°) 


7. Compute tractions: r, = D t] A J = 6 , , K 


1 - d ( 1 + S 3 


Ai 


8. Compute tangent stiffness tensor: 


ta tan 

ij ~ 


Dij - K 




1 + 4^] ^S^a,a,} y<x<Af 


D 




, r* > A or A- f < A 


Further detail regarding the damage model used here can be found in refer- 
ences [2] -[3]. 
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2.2 Cohesive zone model and FEM 


In a finite element model using the CZM approach, the complete material 
description is separated into fracture properties captured by the constitutive 
model of the cohesive surface and the properties of the bulk material, captured 
by the continuum regions. 

There are two conditions to obtain a successful FEM simulation using CZM 
[12]: (a) The cohesive contribution to the global compliance should be small 
enough to avoid the introduction of a fictitious compliance to the model, and 
(b) the element size needs to be less than the cohesive zone length. 

(a) Stiffness of the cohesive zone model 

Different guidelines have been proposed for selecting the stiffness of the inter- 
face. Daudeville et al. [13] calculated the stiffness in terms of the thickness and 
the elastic modulus of the interface. The interface thickness between plies is 
a very small resin-rich region; therefore, the interface stiffness obtained from 
the Daudeville equations is very high. Zou et al. [14] , based on their own expe- 
rience, proposed a value for the interface stiffness between 10 4 and 10' times 
the value of the interfacial strength per unit length. Camanho and Davila [15] 
obtained accurate predictions using a value of 10 6 N/mm 3 . 

The effective clastic properties of the composite depend on the properties of 
both the cohesive surfaces and the bulk constitutive relations. Although the 
compliance of the cohesive surfaces can contribute to the global deformation, 
its only purpose is to simulate fracture. Moreover, the clastic properties of the 
cohesive surfaces are mesh-dependent because the surface relations exhibit an 
inherent length scale that is absent in homogeneous deformations [16]. 

If the cohesive contribution to the compliance is not small enough compared 
to that of the volumetric constitutive relation, a stiff connection between two 
neighboring layers before delamination initiation is not assured. The effect of 
compliance of the interface on the bulk properties of a laminate is illustrated 
in the one- dimensional model shown in Figure 3. The traction continuity con- 
dition requires: 


a = E 3 e = K A (10) 

where a is the traction on the surface, t is the thickness of an adjacent sub- 
laminate, £ = j is the transverse strain, K is the interface stiffness that 
relates the resulting tractions at the interface with the opening displacement 
A, and E% is the through-the-thickness Young’s modulus of the material. For 
a transversely isotropic material E 5 = E 2 . 



The effective strain of the composite is: 



▼ F 


Fig. 3. Influence of the cohesive surface on the deformation. 


( 11 ) 


Since the traction continuity condition requires that a = E e g£ e g, the equivalent 
Young’s modulus E e g can be written as a function of the Young’s modulus of 
the material, the mesh size, and the interface stiffness. Using equations (10) 
and (11), the effective Young’s modulus can be written as: 


E e g — E 3 



( 12 ) 


The effective elastic properties of the composite will not be affected by the 
cohesive surface whenever the inequality E 3 <C Kt is being accomplished, i.e: 


t 

where a is a parameter much larger than 1 (a > 1). However, large values of 
the interface stiffness may cause numerical problems, such as spurious oscilla- 
tions of the tractions [17]. Thus, the interface stiffness should be large enough 
to provide a reasonable stiffness but small enough to avoid numerical problems 
such as spurious oscillations of the tractions in an element. 

The ratio between the value of the Young modulus obtained with equation 
(12) and the Young modulus of the material, as a function of the parameter 
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a is shown in Figure 4. For values of a greater than 50, the loss of stiffness 
due to the presence of the interface is less than 2% (see Fig. 4). 



Fig. 4. Ratio between the equivalent elastic modulus and the Young’s modulus of the 
material, as a function of the parameter a. 

The use of equation (13) is preferable to guidelines presented in previous work 
[13]- [15] because it results from mechanical considerations, and it provides a 
sufficient stiffness (a times the stiffness of the material) while avoiding spurious 
oscillations caused by an excessively stiff interface. The values of the interface 
stiffness obtained with equation (13) and those used by other authors for 
a carbon/epox composite are shown in Table 3. The material’s transverse 
modulus is llkN/mm 2 , its nominal interfacial strength is r° = 45N/mm 2 , and 
a = 50 is selected. 


Table 3. Interface stiffness K proposed by different authors (N/mm 3 ). 


t(mm ) 0.125 

1 

2 

3 

5 

Equation (13) 4.43xl0 6 

5.5x10 s 

2.75x10 s 

1.83x10 s 

1.1x10 s 

Zou et al. [14] 

4.5x10 s 

< K < 

4.5xl0 8 


Camanho and Davila [15] 10 6 

10 6 

10 6 

10 6 

10 6 


Equation (13) gives a range of the interface stiffness between 10 s and 5xl0 6 N/mm 3 
for a sub-laminate thickness between 0.125mm and 5mm. These values are 
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close to the interface stiffness proposed by Camanho and Davila and the val- 
ues obtained with Zou’s guidelines (between 4.5xl0 5 and 4.5xl0 8 N/mm 3 ). 


(b) Length of the cohesive zone 

Under single-mode loading, interface damage initiates at a point where the 
traction reaches the maximum nominal interfacial strength, r°. For mixed- 
mode loading, interface damage onset is predicted by a criterion established 
in terms of the normal and shear tractions. According to Fracture Mechanics, 
cracks propagate when the energy release rate reaches a critical value G c . 
The CZM approach prescribes the interfacial normal and shear tractions that 
resist separation and relative sliding at an interface. The tractions, integrated 
to complete separation, yield the fracture energy release rate, G c . The length 
of the cohesive zone l cz is defined as the distance from the crack tip to the 
point where the maximum cohesive traction is attained (see Figure 5). 


i k 
1 r 

Sublaminate 1 

Sublaminate 2 


rizz: 

. X 

Max. traction 

u 

' Crack tip 




Fig. 5. Length of the cohesive zone. 


Different models have been proposed in the literature to estimate the length 
of the cohesive zone. Irwin [18] estimated the size of the plastic zone ahead 
of a crack in a ductile solid by considering the crack tip zone within which 
the von Mises equivalent stress exceeds the tensile yield stress. Dugdale [4] 
estimated the size of the yield zone ahead of a mode I crack in a thin plate 
of an clastic-perfectly plastic solid by idealizing the plastic region as a narrow 
strip extending ahead of the crack tip that is loaded by the yield traction. 
Barenblatt [5] provided an analogue for ideally brittle materials of the Dugdale 
plastic yield zone analysis. Hui [19] estimated the length of the cohesive zone 
for soft elastic solids, and Falk [12] and Rice [20] estimated the length of the 
cohesive zone as a function of the crack growth velocity. The expressions that 
result from these models for the case of plane stress are presented in Table 4. 
It is assumed that the relation between the critical stress intensity factor K c 
and the critical energy release rate G c can be expressed as = G C E. 

All the models described above to predict the cohesive zone length l cz have 
the form: 
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(14) 


l 


CZ 



where E is the Young modulus of the material, G c is the fracture energy 
release rate, r° is the maximum interfacial strength, and M is a parameter that 
depends on each model. The most commonly used models in the literature are 
Hillerborg’s model [6] and Rice’s model [20]. In these models, the parameter 
M is either close or exactly equal to unity. A summary of the different models 
commonly used in the literature, and the equivalent parameter M for plane 
stress are shown in Table 4. In this paper, Hillcrborg’s model is used. 


For the case of orthotropic materials with transverse isotropy, the value of the 
Young’s modulus in equation (14) is the transverse modulus of the material, 
E 2 . 


Table 4. Length of the cohesive zone and equivalent value of the parameter M. 



^ CZ 

M 

Hui [19] 

2 p G c 
( T 0) 2 

0.21 

Irwin [18] 

1 p Gc 
?r (t 0 ) 2 

0.31 

Dugdale [4], Barenblatt [5] 

7T TP Gc 

0.4 

Rice [20], Falk [12] 

9n 77> G c 
32 ' C/ (r°) 2 

0.88 

Hillcrborg [6] 

p G c 

1 


In order to obtain accurate results using CZM, the tractions in the cohesive 
zone must be represented adequately by the finite element spatial discretiza- 
tion. The number of elements in the cohesive zone is obtained with the equa- 
tion: 

N e = l f (15) 

where l e is the mesh size in the direction of crack propagation. 

When the cohesive zone is discretized by too few elements, the fracture energy 
is not represented accurately and the model does not capture the continuum 
field of a cohesive crack. Therefore, a minimum number of elements, N e , is 
needed in the cohesive zone to get successful FEM results. 

However, the minimum number of elements needed in the cohesive zone is 
not well established: Moes and Bclytschko [21], based on the work of Carpin- 
teri and Cornetti [22], suggest using more than 10 elements. However, Falk 
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et al. [12] used between 2 and 5 elements in their simulations. In the para- 
metric study by Davila and Camanho [23], the minimum element length for 
predicting the delamination in a double cantilever beam (DCB) specimen was 
1mm, which leads, using equation (14), to a length of the cohesive zone of 
3.28mm. Therefore, 3 elements in the cohesive zone were sufficient to predict 
the propagation of delamination in Mode I. 


2. 3 Guidelines for the selection of the parameters of the interface with coarser 
mesh 


One of the drawbacks in the use of cohesive zone models is that very fine 
meshes are needed to assure a reasonable number of elements in the cohesive 
zone. The length of the cohesive zone given by equation (14) is proportional 
to the fracture energy release rate (G c ) and to the inverse of the square of 
the interfacial strength r°. For typical graphite-epoxy or glass-epoxy com- 
posite materials, the length of the cohesive zone is smaller than one or two 
millimeters. Therefore, according to equation (15), the mesh size required in 
order to have more than two elements in the cohesive zone should be smaller 
than half a millimeter. The computational requirements needed to analyze 
a large structure with these mesh sizes may render most practical problems 
intractable. 


Alfano and Crisficld [24] observed that variations of the maximum interfa- 
cial strength do not have a strong influence in the predicted results, but that 
lowering the interfacial strength can improve the convergence rate of the so- 
lution. The result of using a lower interfacial strength is that the length of 
the cohesive zone and the number of elements in the cohesive zone increase. 
Therefore, the representation of the softening response ahead of a crack tip is 
more accurate with a lower interface strength. 


It is possible to develop a strategy to adapt the length of the cohesive zone 
to a given mesh size. The procedure consists of determining the value r° of 
the interfacial strength required for a desired number of elements (Nf ) in the 
cohesive zone. From equations (14) and (15), the required interface strength 
is: 


r 


0 



(16) 


Finally, the interfacial strength is chosen as: 

T = min{T°,T 0 } (17) 


The effect of a reduction of the interfacial strength is to enlarge the cohesive 
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zone, and thus, the model is better suited to capture the softening behaviour 
ahead of the crack tip. Moreover, if equation (13) is used to compute the 
interface stiffness, the interface stiffness will be large enough to assure a stiff 
connection between the two neighboring layers and small enough to avoid 
spurious oscillations. The drawback associated with reducing the interfacial 
strength is that the stress distribution in the regions near the crack tip may 
not be accurate [24], 


3 Simulation of the double cantilever beam specimen 


The influence of mesh size, interface stiffness, and interface strength were in- 
vestigated by analyzing the Mode I test of a double cantilever beam (DCB). 
The DCB specimen was fabricated with a unidirectional T300/977-2 carbon- 
fiber-rcinforced epoxy laminate. The specimen is 150-mm-long, 20.0 nun-wide, 
with two 1.98-mm-thick arms, and an initial crack length of 55mm. The ma- 
terial properties are shown in Table 5 [25]. 

Table 5. Mechanical and interface material properties of T300/977-2. 


An 

E 22 — A33 

G12 — G13 

g 23 

150.0GPa 

ll.OGPa 

6.0GPa 

3.7GPa 

^1 2 = Zb 3 

Z / 23 

Gic 

r° 

'3 

0.25 

0.45 

0.352N/mm 

60MPa 


The FEM model was composed of two layers of four-noded 2D plane strain 
elements connected together with four-node decohesion elements. The decohe- 
sion elements were implemented using a user-written subroutine in the finite 
element code ABAQUS [26]. 

Three sets of simulations were performed. First, a DCB test was simulated 
with different levels of mesh refinement using the material properties shown 
in Table 5 and interfacial stiffness of K=10 6 N/nnn 3 . Then, equations (15) 
and (13) were used to calculate an adjusted interfacial strength and interface 
stiffness. Finally, a set of simulations with a constant mesh size using different 
interface stiffnesses in order to investigate the influence of the stiffness on the 
calculated results. 

Several analyses were carried out for mesh sizes ranging between 0.125mm 
and 5mm. The load-displacement curves obtained for different element sizes 
using the nominal interfacial strength are shown in Figure 6. 
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The results indicate that a mesh size of l e < 0.5mm is necessary to obtain 
converged solutions. The predictions made with coarser meshes significantly 
overpredict the experimental results. 

The length l cz of the cohesive zone for the material given in Table 5 is close 
to 1mm. For a mesh size greater than 0.5mm, fewer than two elements would 
span the cohesive zone, which is not sufficient for an accurate representation 
[22]- [23]. For a mesh size smaller than 0.5mm, more than two elements would 
span the cohesive zone. For a mesh size of 0.25mm, four elements would span 
the cohesive zone. 



Fig. 6. Load- displacement curves obtained for a DCB test with different mesh sizes. 


3. 1 Effect of interface strength 


To verify the effect of interface strength on the computed results, simula- 
tions were performed by specifying the desired number of elements within the 
cohesive zone to be N 0 = 5 and reducing the interface strength according 
to equation (17). The load-displacement curves obtained for several levels of 
mesh refinement are shown in Figure 7. Accurate results are obtained for mesh 
sizes smaller than 2.5nnn. 

A comparison of the load-displacement curves for the DCB specimen calcu- 
lated using the nominal interface strength and the strength obtained from 
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Fig. 7. Load- displacement curves obtained for a DCB test with different mesh sizes 
and keeping constant the number of elements (5) in the cohesive zone. 



Fig. 8. Maximum load obtained in a DCB test for two cases: a) with constant inter- 
facial strength, b) with interfacial strength calculated according to Eq. (17). 
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equation (17) is shown in Figure 8. The maximum load obtained by keep- 
ing the maximum interfacial strength constant increases with the mesh size, 
so the results obtained are mesh dependent, especially for mesh sizes greater 
than 2mm. However, the loads predicted by modifying the interfacial strength 
according to equation (17) are nearly constant for element sizes smaller than 
3mm. In addition, the global deformation and the crack tip position are also 
nearly independent for mesh refinement, as illustrated in Figure 9. 



Fig. 9. Crack tip for different element size. 


3.2 Effect of interface stiffness 


The DCB test was simulated with a mesh size of 2.5mm for various values 
of the interface stiffness in order to investigate the influence of the stiffness 
on the predicted failure load. The results of the simulations are presented in 
Figure 10. 

The load-displacement response curves obtained from simulations using an 
interface stiffness greater than 10 4 N/mm 3 are virtually identical. However, 
smaller values of the interface stiffness have a strong influence on the load- 
displacement curves, since a stiff connection between the two neighboring lay- 
ers is not assured. Moreover, the number of iterations needed for the solution 
when using an interface stiffness smaller than 10 4 N/mm 3 is greater than the 
number of iterations needed for range of the interface stiffness between 10 6 
and 10 10 N/mm 3 . For values of the interface stiffness significantly greater than 
10 lo N/mm 3 , the number of iterations needed for the solution increases (see 
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Figure 11). The stiffness that results from equation (13) is 5.55xl0 5 N/mm 5 , 
which is ideal for good convergence of the solution procedure. 
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Fig. 10. Influence of the value of the interface stiffness on the load- displacement 
curves. 



Fig. 11. Influence of the value of the interface stiffness on the number of iterations. 
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4 Concluding remarks 


An engineering solution for the simulation of delamination using coarse meshes 
was presented. Two new guidelines for the selection of the parameters for the 
constitutive equation used for the simulation of delamination were presented. 

First, a new equation for the selection of the interface stiffness parameter K 
was derived. The new equation is preferable to previous guidelines because 
it results from mechanical considerations rather than from experience. The 
approach provides an adequate stiffness to ensure a sufficiently stiff connec- 
tion between two neighboring layers, while avoiding the possibility of spurious 
oscillations in the solution caused by overly stiff connections. 

Finally, an expression to adjust the maximum interfacial strength used in the 
computations with coarse meshes was presented. It was shown that a minimum 
number of elements within the cohesive zone is necessary for accurate simula- 
tions. By reducing the maximum interfacial strength, the cohesive zone length 
is enlarged and more elements span the cohesive zone. The results obtained 
by reducing the maximum interfacial strength show that accurate results can 
be obtained with a mesh ten times coarser than by using the nominal inter- 
face strength. The drawback in using a reduced interfacial strength value is 
that the stress concentrations in the bulk material near the crack tip are less 
accurate. 
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